home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / SigProc2.0 / Packages / SignalProcessing / Support / TransSupport.m < prev    next >
Encoding:
Text File  |  1992-08-18  |  48.1 KB  |  1,461 lines
  1. (*  :Title:    General Platform for Transforms  *)
  2.  
  3. (*  :Authors:    Brian Evans, James McClellan  *)
  4.  
  5. (*
  6.     :Summary:    Platform for symbolic transforms
  7.         Plots magnitude/phase responses
  8.         Generates pole-zero diagrams and root loci
  9.         Tools for stability analysis based on transform objects.
  10.  *)
  11.  
  12. (*  :Context:    SignalProcessing`Support`TransSupport`  *)
  13.  
  14. (*  :PackageVersion:  2.7    *)
  15.  
  16. (*
  17.     :Copyright:    Copyright 1989-1991 by Brian L. Evans
  18.         Georgia Tech Research Corporation
  19.  
  20.     Permission to use, copy, modify, and distribute this software
  21.     and its documentation for any purpose and without fee is
  22.     hereby granted, provided that the above copyright notice
  23.     appear in all copies and that both that copyright notice and
  24.     this permission notice appear in supporting documentation,
  25.     and that the name of the Georgia Tech Research Corporation,
  26.     Georgia Tech, or Georgia Institute of Technology not be used
  27.     in advertising or publicity pertaining to distribution of the
  28.     software without specific, written prior permission.  Georgia
  29.     Tech makes no representations about the suitability of this
  30.     software for any purpose.  It is provided "as is" without
  31.     express or implied warranty.
  32.  *)
  33.  
  34. (*  :History:    *)
  35.  
  36. (*  :Keywords:    pole-zero diagrams, root locus, frequency response, stability *)
  37.  
  38. (*
  39.     :Source:    {Multidimensional Digital Signal Processing}, 1984,
  40.           by D. Dudgeon and R. Mersereau
  41.         {Digital Signal Processing}, 1975, Oppenheim and Schafer
  42.  *)
  43.  
  44. (*
  45.     :Warning:    The transform objects LTransData and ZTransData must
  46.         be defined before this file is loaded (see "ROC.m").
  47.  *)
  48.  
  49. (*  :Mathematica Version:  1.2 or 2.0  *)
  50.  
  51. (*  :Limitation:  *)
  52.  
  53. (*
  54.     :Discussion:  This file has five sections:
  55.  
  56.         A.  Routines to aid in the finding of symbolic transforms.
  57.         B.  Magnitude and phase plots (for 1-D and 2-D signals)
  58.         C.  Root-locus plots
  59.         D.  Pole-zero diagrams and root maps (for 1-D and
  60.             2-D transforms)
  61.         E.  Stability analysis which rely on transform objects.
  62.  
  63.     Section B and D display the results of transform objects.
  64.     Section E extracts information from transform objects.
  65.  *)
  66.  
  67. (*
  68.     :Functions:    AddT
  69.         ConjT
  70.         ConvolveT
  71.         DerivativeT
  72.         IntegrateT
  73.         LineImpulsemDT
  74.         MagPhasePlot
  75.         MultiDInvTransform
  76.         MultiDLookup
  77.         MultiDTransform
  78.         MultT
  79.         PoleZeroPlot
  80.         ROCPlot
  81.         RootLocus
  82.         ScaleT
  83.         ShadedAnnulus
  84.         Stable
  85.         SubstituteForT
  86.         SubT
  87.         TransformDialogue
  88.         TransformFixUp
  89.  *)
  90.  
  91.  
  92.  
  93. (*  Enhance usage messages  *)
  94.  
  95. If [ TrueQ[ $VersionNumber >= 2.0 ],
  96.      $NewMessage[Apart, "usage"];
  97.      $NewMessage[Definition, "usage"];
  98.      $NewMessage[Simplify, "usage"] ];
  99.  
  100. Apart::usage = Apart::usage ~StringJoin~
  101.    " It is also an option for the inverse z and Laplace transforms \
  102.    directing how to handle taking partial fractions of denominators \
  103.    of 5th degree or less. \
  104.    A value of Rational directs exact roots to be taken whereas All \
  105.    indicates numeric roots."
  106.  
  107. Definition::usage = Definition::usage ~StringJoin~
  108.    " It is also an option for all of the symbolic linear transforms \
  109.    which indicates whether or not to apply the definition of the \
  110.    transform in terms of built-in Mathematica routines if all other
  111.    attempts at the transform have failed."
  112.  
  113. Simplify::usage = Simplify::usage ~StringJoin~
  114.    " It is also an option for z, Laplace, and continuous Fourier transforms, \
  115.    which can take the value of True or False."
  116.  
  117.  
  118. (*  B E G I N     P A C K A G E  *)
  119.  
  120. BeginPackage[ "SignalProcessing`Support`TransSupport`",
  121.           "SignalProcessing`Support`ROC`",
  122.           "SignalProcessing`Support`SigProc`",
  123.           "SignalProcessing`Support`SupCode`" ]
  124.  
  125.  
  126. If [ TrueQ[ $VersionNumber >= 2.0 ],
  127.      Off[ General::spell ];
  128.      Off[ General::spell1 ] ]
  129.  
  130.  
  131. (*  U S A G E     I N F O R M A T I O N  *)
  132.  
  133. AddT::usage =
  134.     "AddT[transq, t1, t2] adds the two transforms t1 and t2 together \
  135.     and determines the new region of convergence. \
  136.     The new transform is returned as a list. \
  137.     AddT[transq, t1, t2, lowerlimit] uses lowerlimit as the lower limit \
  138.     on region of convergence values when combining ROC's (default is 0)."
  139.  
  140. ConjT::usage =
  141.     "ConjT[transq, X(s), s] implements T{ Conj{x(t)} } --> X*(s*), \
  142.     where X(s) is the transform of x(t)."
  143.  
  144. ConvolveT::usage =
  145.     "ConvolveT[transq, convop, t1, t2] implements one dimension of the  \
  146.     transform of a convolution, Convolution[All, All, t][x1, x2], where \
  147.     convop is the convolution operator, t1 is the transform of x1, \
  148.     and t2 is the transform of x2."
  149.  
  150. DerivativeT::usage =
  151.     "DerivativeT[transq, expr, s, m] applies the derivative operator \
  152.     to the transform expr which is a function of s. \
  153.     The derivative is taken m times if and only if expr satisfies transq."
  154.  
  155. DomainScale::usage =
  156.     "DomainScale is an option for MagPhasePlot indicating the \
  157.     scaling of the independent variable axis, Linear or Log, \
  158.     for the magnitude plot."
  159.  
  160. IntegrateT::usage =
  161.     "IntegrateT[transq, t, variable, lower-limit, upper-limit] integrates \
  162.     the transform t (with respect to variable) from lower-limit to \
  163.     upper-limit. \
  164.     The resulting transform is returned as a list."
  165.  
  166. Linear::usage =
  167.     "Linear is a possible value for axis scaling.  See MagPhasePlot."
  168.  
  169. LineImpulsemDT::usage =
  170.     "LineImpulsemDT[transq, expr, s, slist, nleft, fun] applies rules for \
  171.     evaluating multidimensional transforms which involve line impulses, \
  172.     like the z-transform of f[n1,n2] Impulse[n1 - n2] or the Laplace \
  173.     transform of f[t1,t2] Delta[t1 - t2]. \
  174.     In either case, the impulse can be removed by setting n1=n2 or t1=t2, \
  175.     which yields a one-dimensional function. \
  176.     The resulting 1-D transform is then altered by the rule z = z1 z2 \
  177.     or s = s1 + s2, respectively. \
  178.     Here, the argument fun is Times for the z-transform and Plus for \
  179.     the Laplace transform. \
  180.     Note that these rules are only applied if expr is a valid transform, \
  181.     which occurs if transq[expr] evaluates to True."
  182.  
  183. MagRangeScale::usage =
  184.     "MagRangeScale is an option for MagPhasePlot indicating the \
  185.     scaling of the dependent variable axis, Linear or Log, \
  186.     for the magnitude plot. \
  187.     A setting of Null disables the generation of the magnitude plot."
  188.  
  189. MagPhasePlot::usage =
  190.     "MagPhasePlot[freqresp, {w, wmin, wmax }] plots the magnitude \
  191.     and phase response over the specified range of frequencies. \
  192.     It returns a list of two elements:  the magnitude plot and the \
  193.     phase plot. \
  194.     The two-dimensional version has the form \
  195.     MagPhasePlot[freqresp, {w1, wmin1, wmax1}, {w2, wmin2, wmax2}]. \
  196.     The default options are initially biased for continuous-time \
  197.     frequency responses of digital signals. \
  198.     MagnitudePhasePlot is an alias for MagPhasePlot."
  199.  
  200. MultiDInvTransform::usage =
  201.     "MultiDInvTransform[expression, transvar, timevar, options, transq, \
  202.     transform, makeobject, tdefault] finds the multidimensional inverse \
  203.     transform of expression by one call to transform per dimension. \
  204.     Here, transq, transform, and makeobject are function heads."
  205.  
  206. MultiDLookup::usage =
  207.     "MultiDLookup[options, timevars, transformvars] expands any \
  208.     multidimensional transform pairs into a set of one-dimensional \
  209.     transform pairs."
  210.  
  211. MultiDTransform::usage =
  212.     "MultiDTransform[makefun, transform, transtest, expr, expr-vars, \
  213.     def-trans-var, trans-vars] finds the multidimensional transform \
  214.     by calling transform once per dimension. \
  215.     If a call to transform produces an incomplete (invalid) transform, \
  216.     then the multidimensional transform does not exist and this routine is \
  217.     exited. \
  218.     The dimension of the transform is determined from \
  219.     the number of expr-vars (which equals the number of trans-vars). \
  220.     Here, expr is transformed."
  221.  
  222. MultT::usage =
  223.     "MultT[transq, t1, t2] multiplies the transforms t1 and t2 together. \
  224.     MultT[transq, t1, t2, lowerlimit] uses lowerlimit as the lower \
  225.     limit on the region of convergence (ROC) when combining the ROC's of \
  226.     t1 and t2."
  227.  
  228. PhaseRangeScale::usage =
  229.     "PhaseRangeScale is an option for MagPhasePlot indicating the \
  230.     scaling of the dependent variable axis, Linear or Log,
  231.     for the phase plot. \
  232.     A setting of Null disables the generation of the phase plot."
  233.  
  234. PoleZeroPlot::usage =
  235.     "PoleZeroPlot[transform-object] plots the poles and zeros of the \
  236.     transform function as well as the region of convergence. \
  237.     For the z-transform, this function will also plot the unit circle and \
  238.     shade the region of convergence (whenever the ROC does not cover the \
  239.     entire z plane). \
  240.     To plot a transfer function f, use \
  241.     PoleZeroPlot[f, transform-variable, rminus, rplus, zdomainflag]. \
  242.     In the two-dimensional case, the transform-variable, rminus, and \
  243.     rplus are two-element lists. \
  244.     The two-dimensional pole-zero plot has three cases: \
  245.     separable transform, symmetric transform, and root map. \
  246.     The function returns the poles of the transform."
  247.  
  248. Radian::usage =
  249.     "Radian is a possible value for axis scaling.  See MagPhasePlot."
  250.  
  251. ROCPlot::usage =
  252.     "ROCPlot[{rm, rp}] will plot the region of convergence as a \
  253.     shaded annular region. \
  254.     Possible options for ROCPlot are the same as those for the \
  255.     Show (or Plot) command."
  256.  
  257. RootLocus::usage =
  258.     "RootLocus[f(z), z, {freeparam, start, end, step}] plots \
  259.     the root locus of f(z) with respect to freeparam over the \
  260.     range of start to end at increments of step. \
  261.     In order to connect the points, the function f(z) will be factored \
  262.     symbolically with respect to z. \
  263.     All closed-form factors of reasonable size will be graphed as \
  264.     parametric plots. \
  265.     RootLocus will generate three plots:  the pole zero plot for \
  266.     freeparam = start, the root locus, and the pole-zero plot for \
  267.     freeparam = end. \
  268.     The function returns a list of all three plots. \
  269.     RootLocus supports the same options as does Show. \
  270.     RootLocusPlot is an alias of RootLocus."
  271.  
  272. ScaleT::usage =
  273.     "ScaleT[transq, t, c] multiplies the transform t by c while leaving \
  274.     the region of convergence unaltered. \
  275.     The resulting transform is returned as a list."
  276.  
  277. ShadedAnnulus::usage =
  278.     "ShadedAnnulus[rm, rp] and ShadedAnnulus[rm, rp, angulartilt] \
  279.     will create a 2-D graphics object that is an annulus \
  280.     (rm < radius < rp) filled with uniformly spaced horizontal \
  281.     lines rotated by angulartilt degrees."
  282.  
  283. Stable::usage =
  284.     "Stable[transexpr] returns True if the transform transexpr \
  285.     represents a stable time-domain signal."
  286.  
  287. SubT::usage =
  288.     "SubT[transq, t1, t2] subtracts the transforms t2 from t1 and \
  289.     determines the new region of convergence. \
  290.     The resulting transform is returned as a list. \
  291.     SubT[transq, t1, t2, lowerlimit] uses lowerlimit as the lower limit \
  292.     on region of convergence values when combining ROC's (default is 0)."
  293.  
  294. SubstituteForT::usage =
  295.     "SubstituteForT[transq, t, s, news] substitutes news at every \
  296.     occurrence of s in the transform t. \
  297.     The resulting transform is returned as a list."
  298.  
  299. TransformDialogue::usage =
  300.     "TransformDialogue[ strategy, options, arg1, arg2, ... ] is the \
  301.     routine that displays textual dialogue about the strategy being used \
  302.     by a particular transform."
  303.  
  304. TransformFixUp::usage =
  305.     "TransformFixUp[time_var, transform_var, options, \
  306.     transform_head, triplet_flag, transform_name, lower, upper] \
  307.     attempts to resolve incomplete transforms according to the \
  308.     transform pairs provided by the user via the TransformLookup option. \
  309.     If the transform is a triplet of information (funct, Rminus, Rplus), \
  310.     e.g. the forward z- and Laplace transforms, then triplet_flag is True \
  311.     and the lower and upper arguments are used for the second and \
  312.     third fields of the triplet."
  313.  
  314. TransformLookup::usage =
  315.     "TransformLookup is an option for the transform rule bases. \
  316.     It allows the user to specify a list of transform pairs to augment \
  317.     those that already exist in the transform rule bases. \
  318.     For example, ZTransform[ Shift[L, n][x[n]], n, z, \
  319.     TransformLookup -> { x[n] :> X[z] } ]. \
  320.     To specify a region of convergence, use \
  321.     TransformLookup -> { x[n] :> { X[z], rm, rp } }. \
  322.     This example allows the user to work x[n] without every defining it. \
  323.     The multidimensional case is more tricky: \
  324.     ZTransform[ x[n1, n2], n, z, TransformLookup -> \
  325.     { x[n1, n2] :> X[z1, n2], X[z1, n2] :> X[z1, z2] } ]. \
  326.     A region of convergence can be specified for each dimension. \
  327.     For the z- and Laplace transform rule bases, one can also specify \
  328.     a region of convergence associated with the transform pairs."
  329.  
  330. (*  E N D     U S A G E     I N F O R M A T I O N  *)
  331.  
  332. If [ TrueQ[ $VersionNumber >= 2.0 ],
  333.      Terms::usage = "An"; Protect[ Terms ],
  334.      Terms::usage = Terms::usage ~StringJoin~ " It is also an" ]
  335.  
  336. Terms::usage = Terms::usage ~StringJoin~
  337.     " option for the inverse z- and Laplace transform rule bases, \
  338.     InvZTransform and InvLaPlace, which specifies that the inverse \
  339.     will be approximated by a series expansion if all other attempts \
  340.     at the inverse have failed.";
  341.  
  342.  
  343. Begin["`Private`"]
  344.  
  345.  
  346. (*  C O N S T A N T S  *)
  347.  
  348. axesdefaultvalue = If [ TrueQ[$VersionNumber >= 2.0], True, Automatic ]
  349.  
  350.  
  351. (*  M E S S A G E S  *)
  352.  
  353. MagPhasePlot::unresolved =
  354.     "These symbols were assigned a value of 1:  ``."
  355. MagPhasePlot::nodeltas =
  356.     "Delta functions are not shown on the graph"
  357. MultiDInvTransform::notenough =
  358.     "Not enough non-transform (time) variables specified."
  359. MultiDTransform::notenough =
  360.     "Not enough transform domain variables specified."
  361. PoleZeroPlot::invalidROC =
  362.     "The poles `` are inside the specified region of convergence {``, ``}."
  363. PoleZeroPlot::multidimensions =
  364.     "Can not produce a pole-zero plot for the multidimensional transform."
  365. PoleZeroPlot::notrational =
  366.     "Transform is not a rational polynomial."
  367. PoleZeroPlot::noplot =
  368.     "A pole-zero plot cannot be generated."
  369. PoleZeroPlot::posROC =
  370.     "Possible Regions of Convergence are: "
  371. RootLocus::constant =
  372.     "The function `` does not depend on the variable ``."
  373. RootLocus::notnum =
  374.     "The min, max, or step value is not a number: ``."
  375. Transform::badvar =
  376.     "The `` variable given to `` is not a symbol or list of symbols: ``"
  377. Transform::incomplete =
  378.     "The rule base could not compute the `` of `` with respect to ``."
  379. Transform::integrate =
  380.     "Integrating the function `` with respect to ``...."
  381. Transform::novariables =
  382.     "No `` variables specified for transform.  Possible variables are ``."
  383. Transform::notenough =
  384.     "Not enough `` variables specified for the transform."
  385. Transform::twosided =
  386.     "The two-sided transform was applied to ``."
  387. TransformLookup::nomatch
  388. TransformLookup::notrule =
  389.     "The expression `` passed in the TransformLookup option of the `` \
  390.     object is not a rule."
  391. TransformFixUp::notriplet =
  392.     "The rule `` has been changed to `` because `` does not track \
  393.      a region of convergence."
  394.  
  395.  
  396. (*  A.  B A S I C     R O U T I N E S     F O R     T R A N S F O R M S  *)
  397.  
  398. (*  AddT  *)
  399. AddT[transq_, t1_, t2_, lowerlimit_:0] :=
  400.     Transform[ TheFunction[t1] + TheFunction[t2],
  401.            FindRMinus[ GetRMinus[t1], GetRMinus[t2], lowerlimit ],
  402.            FindRPlus[ GetRPlus[t1], GetRPlus[t2], lowerlimit ] ] /;
  403.     transq[t1] && transq[t2]
  404.  
  405. Format[ AddT[transq_, t1_, t2_, lowerlimit_:0] ] := t1 + t2
  406.  
  407. (*  ConjT  *)
  408. ConjT[transq_, t_, s_] :=
  409.     Transform[ Conjugate[TheFunction[t]] /. s -> Conjugate[s],
  410.            GetRMinus[t], GetRPlus[t] ] /;
  411.     transq[t]
  412.  
  413. Format[ ConjT[transq_, t_, s_] ] := 
  414.     SequenceForm[ ToString[StringForm["{ Conj{``} }", t]],
  415.               Subscript[s -> news] ]
  416.  
  417. (*  ConvolveT  *)
  418. ConvolveT[transq_, convop_, t1_, t2_] :=
  419.     Block [    {result},
  420.         result = MultT[transq, t1, t2];
  421.         Transform[ convop [t1, t2], GetRMinus[t], GetRPlus[t] ] ] /;
  422.     transq[t1] && transq[t2]
  423.  
  424. Format[ ConvolveT[transq_, convop_, t1_, t2_] ] :=
  425.     ToString[StringForm["`` ** ``", t1, t2]]
  426.  
  427. (*  DerivativeT  *)
  428. DerivativeT[transq_, x_, s_, m_] :=
  429.     Transform[ D[TheFunction[x], {s, m}], GetRMinus[x], GetRPlus[x] ] /;
  430.     transq[x]
  431.  
  432. Format[ DerivativeT[transq_, x_, s_, m_] ] :=
  433.     SequenceForm[ ColumnForm[{"D" Superscript[m],
  434.                   "  " ~StringJoin~ ToString[s]}],
  435.               { x } ]
  436.  
  437. (*  IntegrateT  *)
  438. IntegrateT[transq_, t_, var_, lower_, upper_] :=
  439.     Transform[ Integrate[ TheFunction[t], {var, lower, upper} ],
  440.            GetRMinus[t], GetRPlus[t] ] /;
  441.     transq[t]
  442.  
  443. Format[ IntegrateT[transq_, t_, var_, l_, u_] ] :=
  444.     ToString[StringForm["Integrate[``, {``, ``, ``}]", t, var, l, u]]
  445.  
  446. (*  LineImpulsemDT  *)
  447. LineImpulsemDT[transq_, t_, s_, slist_, nleft_, fun_] :=
  448.     Block [    {},
  449.         Needs[ "SignalProcessing`Support`MDSupport`" ];
  450.         lineImpulsemDT[transq, t, s, slist, nleft, fun] ] /;
  451.     transq[t]
  452.  
  453. Format[ SignalProcessing`ROCinfo[n_, rm_, rp_] ] :=
  454.     Subscript[ "ROC", {rm, rp} ]
  455.  
  456. (*  MultiDInvTransform  *)
  457. MultiDInvTransform[e_, s_, t_, op_, transq_, trans_, makeobject_, tdefault_] :=
  458.     Block [ {},
  459.         Needs[ "SignalProcessing`Support`MDSupport`" ];
  460.         multiDInvTransform[e, s, t, op, transq,
  461.                    trans, makeobject, tdefault] ]
  462.  
  463. (*  MultiDLookup  *)
  464. MultiDLookup[op_, nlist_, zlist_] :=
  465.     Block [ {},
  466.         Needs[ "SignalProcessing`Support`MDSupport`" ];
  467.         multiDLookup[op, nlist, zlist] ]
  468.  
  469. (*  MultiDTransform  *)
  470. MultiDTransform[makeobjectfun_, transform_, transq_, e_, time_, rest___] :=
  471.     Block [ {},
  472.         Needs[ "SignalProcessing`Support`MDSupport`" ];
  473.         multiDTransform[makeobjectfun, transform,
  474.                 transq, e, time, rest] ]
  475.  
  476. (*  MultT  *)
  477. MultT[transq_, t1_, t2_, lowerlimit_:0] :=
  478.     Transform[ TheFunction[t1] TheFunction[t2],
  479.            FindRMinus[ GetRMinus[t1], GetRMinus[t2], lowerlimit ],
  480.            FindRPlus[ GetRPlus[t1], GetRPlus[t2], lowerlimit ] ] /;
  481.     transq[t1] && transq[t2]
  482.  
  483. Format[MultT[transq_, t1_, t2_, lowerlimit_:0]] := t1 t2
  484.  
  485. (*  ScaleT  *)
  486. ScaleT[transq_, t_, c_] :=
  487.     Transform[ c TheFunction[t], GetRMinus[t], GetRPlus[t] ] /;
  488.     transq[t]
  489.  
  490. ScaleT[transq_, ScaleT[transq_, t_, c_], d_] := ScaleT[transq, t, c d]
  491.  
  492. Format[ScaleT[transq_, t_, c_]] := c t
  493.  
  494. (*  SubstituteForT  *)
  495. SubstituteForT[transq_, t_, s_, news_] :=
  496.     Transform[ TheFunction[t] /. s -> news, GetRMinus[t], GetRPlus[t] ] /;
  497.     transq[t]
  498.  
  499. Format[SubstituteForT[transq_, t_, s_, news_]] :=
  500.     SequenceForm[ {t}, Subscript[s -> news] ]
  501.  
  502. (*  SubT  *)
  503. SubT[transq_, t1_, t2_, lowerlimit_:0] :=
  504.     Transform[ TheFunction[t1] - TheFunction[t2],
  505.            FindRMinus[ GetRMinus[t1], GetRMinus[t2], lowerlimit ],
  506.            FindRPlus[ GetRPlus[t1], GetRPlus[t2], lowerlimit ] ] /;
  507.     transq[t1] && transq[t2]
  508.  
  509. Format[SubT[transq_, t1_, t2_, lowerlimit_:0]] := t1 - t2
  510.  
  511. (*  Transform  *)
  512. Transform/: TheFunction[ Transform[fun_, rest___] ] := fun
  513.  
  514. (*  TransformDialogue  *)
  515. TransformDialogue[ Definition, options_, operator_, oldexpr_, trans_ ] :=
  516.     Block [ {},
  517.         If [ InformUserQ[options],
  518.              Print[ "( after using ", operator,
  519.                 "  to apply the definition to" ];
  520.              Print[ "  ", oldexpr ];
  521.              Print[ "  to find its transform"];
  522.              Print[ "  ", trans, " )" ] ];
  523.         trans ]
  524.  
  525.  
  526. (*  TransformFixUp: resolves functions according to user's transform pairs *)
  527. TransformFixUp[ n_, z_, {}, head_, triplet_, name_, low_, up_ ] := {}
  528.  
  529. makefixuprules[n_, triplet_, name_] :=
  530.       {    (a_ -> b_) :> ToCollection[{}] /; FreeQ[a, n],
  531.     (a_ :> b_) :> ToCollection[{}] /; FreeQ[a, n],
  532.     If [ TrueQ[triplet],
  533.          (a_ -> {b_, rm_, rp_}) :> (a -> Transform[b, rm, rp]),
  534.          (a_ -> {b_, rm_, rp_}) :>
  535.         MyMessage[ TransformFixUp::notriplet, a -> b,
  536.                a -> {b, rm, rp}, a -> b, name ] ],
  537.     If [ TrueQ[triplet],
  538.          (a_ :> {b_, rm_, rp_}) :> (a :> Transform[b, rm, rp]),
  539.          (a_ :> {b_, rm_, rp_}) :>
  540.         MyMessage[ TransformFixUp::notriplet, a :> b,
  541.                a :> {b, rm, rp}, a :> b, name ] ],
  542.     x_ :> x }
  543.  
  544. TransformFixUp[ n_, z_, options_, head_, triplet_, name_, low_, up_ ] :=
  545.     Block [    {fixrules, nolookup, op, rules = {}, validop},
  546.         op = Replace[TransformLookup, options];
  547.         nolookup = SameQ[op, {}] || SameQ[op, TransformLookup];
  548.         If [ ! nolookup,
  549.              fixrules = makefixuprules[n, triplet, name];
  550.              validop = Map[ Replace[#1, fixrules]&, ToList[op] ];
  551.              If [ ! EmptyQ[validop],
  552.               rules = makepostrules[ validop, n, z, head,
  553.                          triplet, name, low, up ] ] ];
  554.  
  555.         rules ]
  556.  
  557. TransformFixUp[ trans_, n_, z_, options_, head_, triplet_, name_, low_, up_ ] :=
  558.     Block [    {rules},
  559.  
  560.         rules = TransformFixUp[ n, z, options, head,
  561.                     triplet, name, low, up ];
  562.         If [ EmptyQ[rules],
  563.              trans,
  564.              ReplaceRepeated[trans, rules] ] ]
  565.  
  566.  
  567. SetAttributes[makepostrules, Listable]
  568.  
  569.     (* rules for transforms returning triplets (forward z- and Laplace) *)
  570. makepostrules[ a_ -> b_Transform, n_, z_, h_, True, transname_, low_, up_ ] :=
  571.     ( h[ a, n, z, __ ] :> b )
  572.  
  573. makepostrules[ a_ -> b_, n_, z_, h_, True, transname_, low_, up_ ] :=
  574.     ( h[ a, n, z, __ ] :> Transform[b, low, up] )
  575.  
  576. makepostrules[ a_ :> b_Transform, n_, z_, h_, True, transname_, low_, up_ ] :=
  577.     ( h[ a, n, z, __ ] :> b )
  578.  
  579. makepostrules[ a_ :> b_, n_, z_, h_, True, transname_, low_, up_ ] :=
  580.     ( h[ a, n, z, __ ] :> Transform[b, low, up] )
  581.  
  582.     (* rules for transforms not returning triplets *)
  583. makepostrules[ a_ -> b_, n_, z_, h_, False, transname_, low_, up_ ] :=
  584.     ( h[ a, n, z, __ ] :> b )
  585.  
  586. makepostrules[ a_ :> b_, n_, z_, h_, False, transname_, low_, up_ ] :=
  587.     ( h[ a, n, z, __ ] :> b )
  588.  
  589.     (* invalid rule encountered *)
  590. makepostrules[ x_, n_, z_, h_, triplet_, transname_, low_, up_ ] :=
  591.     Message[ TransformLookup::notrule, x, transname ]
  592.  
  593.  
  594.  
  595.  
  596. (*  B.  F R E Q U E N C Y     R E S P O N S E  *)
  597.  
  598. (*  getlocation -- Delta functions, a special case of continuous plotting *)
  599.  
  600. getlocation[ c_. Delta[a_. t_ + b_.], t_ ] := N[ - b / a ]
  601.  
  602.  
  603. (*  MagPhasePlot  *)
  604.  
  605. MagPhasePlot/: Options[MagPhasePlot] :=
  606.     { DisplayFunction :> $DisplayFunction,
  607.       Domain -> Continuous, DomainScale -> Linear,
  608.       MagRangeScale -> Linear, PhaseRangeScale -> Degree,
  609.       PlotRange -> All }
  610.  
  611. badoptions[ MagPhasePlot ] =
  612.     { Domain, DomainScale, MagRangeScale, PhaseRangeScale }
  613.  
  614. (*  Supporting local functions  *)
  615.  
  616. fmag[w0_Real, omega_, fresp_, logflag_] :=
  617.     Block [    {value},
  618.         value = Abs[ GetValue[fresp, omega, w0] ];
  619.         If [ logflag, 20 Log[10, value], value ] ]
  620.  
  621. fmag[w01_Real, omega1_, w02_Real, omega2_, fresp_, logflag_] :=
  622.     Block [    {value},
  623.         value = Abs[ GetValue[fresp, {omega1, omega2}, {w01, w02}] ];
  624.         If [ logflag, 20 Log[10, value], value ] ]
  625.  
  626. fphase[w0_Real, omega_, fresp_, indegrees_] :=
  627.     Block [    {value},
  628.         value = GetValue[fresp, omega, w0];
  629.         If [ ZeroQ[value], Return[0.] ];
  630.         value = N [ Arg[value] ];
  631.         Chop[ If [ indegrees, value / Degree, value ] ] ]
  632.  
  633. fphase[w01_Real, omega1_, w02_Real, omega2_, fresp_, indegrees_] :=
  634.     Block [    {value},
  635.         value = GetValue[fresp, {omega1, omega2}, {w01, w02}];
  636.         If [ ZeroQ[value], Return[0.] ];
  637.         value = N [ Arg[value] ];
  638.         Chop[ If [ indegrees, value / Degree, value ] ] ]
  639.  
  640. waxis[w_Real, logflag_] := If [ logflag, Log[10, w], w ]
  641.  
  642. loglabel[w_] := ToString[StringForm["````", Subscripted[Log[10]], w]]
  643.  
  644.  
  645. (*  Supporting rule bases for MagPhase Plot  *)
  646. (*  Interprets signal processing expressions so that they can be plotted  *)
  647.  
  648. ScaleAxisRules = {
  649.     ScaleAxis[sc_, w_][x_] :>
  650.         Block [ {i, result},
  651.             result = x /. w -> sc w;
  652.             For [ i = 1, i < sc, i++,
  653.                   result += x /. w -> (sc w + i 2 Pi) ];
  654.             result ]
  655. }
  656.  
  657. PeriodicRules = {
  658.  
  659.     Periodic[k_, w_Symbol][freqresp_] :>
  660.         Block [ {fresp},
  661.             fresp = TheFunction[freqresp /. ScaleAxisRules];
  662.             fresp + ( fresp /. w -> w + k ) +
  663.                 ( fresp /. w -> w - k ) ],
  664.  
  665.     Periodic[{k1_, k2_}, {w1_Symbol, w2_Symbol}][freqresp_] :>
  666.         Block [    {dim1, dim2, fresp, rules},
  667.             fresp = TheFunction[freqresp /. ScaleAxisRules];
  668.             rules = { w1 -> w1 + dim1 k1,
  669.                   w2 -> w2 + dim2 k2 };
  670.             Apply[Plus,
  671.                   Table[fresp /. rules,
  672.                     {dim1, -1, 1}, {dim2, -1, 1}]] ]
  673.  
  674. }
  675.  
  676.  
  677. (*  Magnitude and phase plots for 2-D signals (must be defined before 1-D)  *)
  678. MagPhasePlot[freqresp_, {w1_Symbol, wmin1_, wmax1_},
  679.                 {w2_Symbol, wmin2_, wmax2_}, op___] :=
  680.     Block [    {},
  681.         Needs[ "SignalProcessing`Support`MDPlotting`" ];
  682.         magPhasePlot2D[ freqresp,
  683.                 {w1, wmin1, wmax1},
  684.                 {w2, wmin2, wmax2},
  685.                 op] ]
  686.  
  687. (*  Magnitude and phase plots for 1-D signals  *)
  688. MagPhasePlot[ freqresp_, {w_Symbol, wminimum_, wmaximum_}, op___ ] :=
  689.     Block [    {arglist, bogus, degrees, deltafuns, deltaplot, disp, fresp,
  690.          ftemp, i, list, logdomain, logrange, magfun,
  691.          magplot = NullPlot, max, maxdec, min, mindec, omega,
  692.          omitplot, options, period, phasefun, phaseplot = NullPlot,
  693.          plotlabel, plotoptions, plotrange, points, ticks,
  694.          varlist, varlist2, wmin, wmax, wstr, xlabel},
  695.  
  696.         (*  Check for errors in arguments  *)
  697.  
  698.         If [ N [ wminimum > wmaximum ],
  699.              Message[Plot::limits, {w, wminimum, wmaximum}];
  700.              Return[Null] ];
  701.  
  702.         (*  Set up for plotting  *)
  703.  
  704.         Off[Power::infy, Infinity::indt];
  705.         options = ToList[op] ~Join~ Options[MagPhasePlot];
  706.         wstr = StripPackage[w];
  707.         plotoptions = RemoveOptions[options, badoptions[MagPhasePlot]];
  708.  
  709.         (*  Separate functions into delta functions and other  *)
  710.         (*  functions keeping in mind the interval of interest *)
  711.         {deltafuns, fresp} =
  712.             SignalCleanup[freqresp, w, wminimum, wmaximum];
  713.  
  714.         (*  If frequency response is zero, plot deltas and exit  *)
  715.  
  716.         If [ ZeroQ[fresp],
  717.              arglist = Join[ plotoptions, 
  718.                      { Axes -> axesdefaultvalue,
  719.                        AxesLabel -> { wstr, " " },
  720.                        PlotLabel -> "Magnitude Response",
  721.                        PlotRange -> All,
  722.                        Ticks -> Automatic } ];
  723.              PrependTo[ arglist,
  724.                 DeltaPlot[ deltafuns, w, N[wminimum],
  725.                        N[wmaximum], 0.0, 1.0 ] ];
  726.              magplot = Apply[ Show, arglist ];
  727.              Return[ { magplot, phaseplot } ] ];
  728.  
  729.         (*  Rename the frequency variable  *)
  730.  
  731.         fresp = fresp /. ( w -> omega );
  732.  
  733.         (*  Find valueless symbols other than frequency variable *)
  734.                 (*  The first parameter in Summation is a local symbol   *)
  735.  
  736.         ftemp = fresp /. 
  737.                           ( Summation[n_, l_, u_, inc_][expr_] :>
  738.                 bogus[l, u, inc, GetVariables[expr /. n -> l]] );
  739.         varlist = Select[GetVariables[ftemp], (! SameQ[#1, omega])&];
  740.         If [ ! EmptyQ[varlist],
  741.              Message[ MagPhasePlot::unresolved, Sort[ varlist ] ];
  742.              fresp = fresp /. Map[(#1 -> 1)&, varlist] ];
  743.  
  744.         (*  Set up for plotting  *)
  745.  
  746.         logrange = SameQ[Replace[MagRangeScale, options], Log];
  747.         plotlabel = If [ logrange,
  748.                  "Magnitude Response (dB)",
  749.                  "Magnitude Response" ];
  750.  
  751.         logdomain = SameQ[Replace[DomainScale, options], Log];
  752.         If [ ZeroQ[fresp], logdomain = False ];
  753.         xlabel = If [ logdomain, loglabel[wstr], wstr ];
  754.         If [ logdomain,
  755.              Off[General::dby0];
  756.              mindec = Floor[ N[ Log[10, wminimum] ] ];
  757.              maxdec = Ceiling[ N[ Log[10, wmaximum] ] ];
  758.              wmin = 10^mindec;
  759.              wmax = 10^maxdec;
  760.              ticks = {};
  761.              For [ i = mindec, i < maxdec, i++,
  762.                ticks = ticks ~Join~
  763.                    {i, i + 0.3, i + 0.7} ];
  764.              ticks = { ticks, Automatic };
  765.              On[General::dby0],
  766.  
  767.              wmin = wminimum;
  768.              wmax = wmaximum;
  769.              ticks = Automatic ];
  770.  
  771.         (*  Magnitude plot with Delta functions (if any)  *)
  772.  
  773.         omitplot = SameQ[Replace[MagRangeScale, options], Null];
  774.         If [ ! omitplot,
  775.              arglist = Join[ { { waxis[w, logdomain],
  776.                          fmag[w, omega, fresp, logrange] },
  777.                        { w, wmin, wmax },
  778.                        DisplayFunction -> Identity },
  779.                      plotoptions ];
  780.              magplot = Apply[ ParametricPlot, arglist ];
  781.  
  782.              (*  determine the range of the magnitude plot  *)
  783.  
  784.              plotrange = Replace[PlotRange, options];
  785.              list = Map[ Second, magplot[[1]] [[1]] [[1]] [[1]] ];
  786.              list = Select[list, NumberQ];    (* remove infinities *)
  787.              max = Max[list];
  788.              min = Min[list];
  789.  
  790.              If [ logrange,
  791.                   max = 20 Ceiling[max / 20];
  792.                   min = 20 Floor[min / 20];
  793.                   If [ ! SameQ[deltafuns, Null] && max < 20,
  794.                    max = 20 ];
  795.                   If [ ( max - min ) > 100,
  796.                    min = max - 100 ];
  797.               If [ SameQ[plotrange, All],
  798.                    plotrange = {min, max} ] ];
  799.  
  800.              If [ ! logrange && max == 0, max = 1 ];
  801.  
  802.              (*  complete magnitude plot --  add in delta functions  *)
  803.  
  804.              magplot =
  805.                Apply[ Show,
  806.                   Join[ { DeltaPlot[deltafuns, omega, wmin, wmax,
  807.                         magplot, min, max] },
  808.                     plotoptions,
  809.                     { PlotLabel -> plotlabel,
  810.                       PlotRange -> plotrange,
  811.                       Ticks -> ticks,
  812.                       AxesLabel -> { xlabel, " " } } ] ] ];
  813.  
  814.         (*  Phase plot  *)
  815.  
  816.         omitplot = SameQ[Replace[PhaseRangeScale, options], Null];
  817.         If [ ! omitplot,
  818.              degrees = SameQ[Replace[PhaseRangeScale, options], Degree];
  819.              plotlabel = If [ degrees,
  820.                       "Phase Response (degrees)",
  821.                       "Phase Response (radians)" ];
  822.              arglist = Join[ { { waxis[w, logdomain],
  823.                      fphase[w, omega, fresp, degrees] },
  824.                        { w, wmin, wmax } },
  825.                     { PlotRange -> All },
  826.                      plotoptions,
  827.                      { PlotLabel -> plotlabel,
  828.                        Ticks -> ticks,
  829.                        AxesLabel -> { xlabel, " " } } ];
  830.              phaseplot = Apply[ ParametricPlot, arglist ] ];
  831.  
  832.         (*  Clean up and return the two plots as graphics objects  *)
  833.         On[Power::infy, Infinity::indt];
  834.         { magplot, phaseplot } ]
  835.  
  836.  
  837. (*  C.  R O O T - L O C U S     P L O T T I N G  *)
  838.  
  839. Options[RootLocus] :=
  840.   { Axes -> axesdefaultvalue, AxesLabel -> { "Re", "" },
  841.     DisplayFunction :> $DisplayFunction, PlotRange -> All }
  842.  
  843. initPlot[ startpoly_, endpoly_, z_, text_, multtext_ ] :=
  844.     Block [    {coords, endplot, rootlist, startplot},
  845.          rootlist = GetRootList[ startpoly, z, N ];
  846.         coords = ComplexTo2DCoordList[rootlist];
  847.          startplot = PointwisePlot[ coords, text, multtext ];
  848.          rootlist = GetRootList[ endpoly, z, N ];
  849.         coords = ComplexTo2DCoordList[rootlist];
  850.          endplot = PointwisePlot[ coords, text, multtext ];
  851.         {startplot, endplot} ]
  852.  
  853. localplot[ label_, oplist_, plots___ ] :=
  854.     Apply[ Show, Join[ { plots }, { PlotLabel -> label }, oplist ] ]
  855.  
  856. oneBranch[ rule_, ktemp_, 0, l_. Pi, kstep_, z_, psize_ ] :=
  857.     Block [ {fun, simprule},
  858.         simprule = simplifyrule[ktemp];
  859.         fun = Replace[z, rule] /. simprule;
  860.         fun = ( fun /. ktemp -> (Pi ktemp) ) /. simprule;
  861.         ParametricPlot[ reImPair[fun],
  862.                 {ktemp, 0, l},
  863.                 PlotStyle -> Thickness[psize/3],
  864.                 DisplayFunction -> Identity ] ]
  865.  
  866. oneBranch[ rule_, ktemp_, kmin_, kmax_, kstep_, z_, psize_ ] :=
  867.     Block [ {fun},
  868.         fun = Replace[z, rule] /. simplifyrule[ktemp];
  869.         ParametricPlot[ reImPair[fun],
  870.                 {ktemp, kmin, kmax},
  871.                 PlotStyle -> Thickness[psize/3],
  872.                 DisplayFunction -> Identity ] ]
  873.  
  874. oneLocus[ poly_, ktemp_, kmin_, kmax_, kstep_, z_, psize_ ] :=
  875.     Block [ {coords, roots},
  876.         roots = Table[ ToCollection[ GetRootList[poly, z, N] ],
  877.                    {ktemp, kmin, kmax, kstep} ];
  878.         coords = Map[Point, ComplexTo2DCoordList[roots]];
  879.         Show[ Prepend[ Map[ oneBranch[#, ktemp, kmin,
  880.                           kmax, kstep, z, psize]&,
  881.                     Select[ Solve[poly == 0, z], tractableQ ] ],
  882.                    Graphics[ Prepend[coords, PointSize[psize]] ] ],
  883.               DisplayFunction -> Identity ] ]
  884.  
  885. polyorder[poly_, z_] := Exponent[poly, z] + Exponent[poly /. z -> z^-1, z]
  886.  
  887. reImPair[f_] := { Re[f], Im[f] }
  888.  
  889. simplifyrule[ ktemp_ ] := 
  890.     Block [    {},
  891.         Exp[ Complex[0, b_] Pi c_. ktemp ] :>
  892.           Exp[ I Mod[b c, 2] Pi ktemp ] /;
  893.           RationalQ[b c] || RealQ[b c] ]
  894.  
  895. (*  Tractable means a closed-form solution of less than 30 terms  *)
  896. tractableQ[ ToRules[x_] ] := False
  897. tractableQ[ {z_ -> a_} ] := ( Count[a, b_, Infinity] < 30 )
  898. tractableQ[ x_ ] := False
  899.  
  900. RootLocus[ poly_, z_Symbol, {k_Symbol, start_, end_, step_:1}, op___ ] :=
  901.     MyMessage[ RootLocus::constant,
  902.            { NullPlot, NullPlot, NullPlot }, poly, z ] /;
  903.     FreeQ[poly, z]
  904.  
  905. RootLocus[ poly_, z_Symbol, {k_Symbol, start_, end_, step_:1}, op___ ] :=
  906.     Block [    {Epplot = NullPlot, Ezplot = NullPlot, Szplot = NullPlot,
  907.          Spplot = NullPlot, denom, expdenom, expnumer, factor,
  908.          kmax, kmin, kstep, ktemp, newpoly, numer, options,
  909.          order, pointthickness, pplot = NullPlot, zplot = NullPlot},
  910.  
  911.         (*  Check numeric arguments  *)
  912.         If [ ! Apply[And, Map[RealValuedQ, N[{start, end, step}]]],
  913.              Message[ RootLocus::notnum, {start, end, step} ];
  914.              Return[ NullPlot ] ];
  915.  
  916.         (*  Disable compilation warnings  *)
  917.         If [ TrueQ[ $VersionNumber >= 2.0 ],
  918.              Off[ ParametricPlot::ppcom ] ];
  919.  
  920.         (*  Determine the numeric range for the free param k  *)
  921.         options = Flatten[ToList[op]] ~Join~ Options[RootLocus];
  922.         If [ N[start < end],
  923.              kmin = start; kmax = end,
  924.              kmin = end; kmax = start ];
  925.         kstep = If [ Positive[step], step, -step ];
  926.  
  927.         (*  Put the function over a common denominator *)
  928.         newpoly = Together[poly] /. k -> ktemp;
  929.  
  930.         (*  Get numerator and denominator--  expand them for speed  *)
  931.         (*  note: expand also hides errors in Solve under Mma 1.2   *)
  932.         expnumer = Expand[Numerator[newpoly]];
  933.         expdenom = Expand[Denominator[newpoly]];
  934.         If [ TrueQ[ $VersionNumber >= 2.0 ],
  935.              numer = Factor[expnumer]; denom = Factor[expdenom],
  936.              numer = expnumer; denom = expdenom ];
  937.  
  938.         (*  Compute the thickness of the points in the locus  *)
  939.         order = 1;
  940.         If [ MixedPolynomialQ[expnumer, z],
  941.              order += polyorder[expnumer, z] ];
  942.         If [ MixedPolynomialQ[expdenom, z],
  943.              order += polyorder[expdenom, z] ];
  944.         factor = N[ Log[10, order (kmax - kmin) / kstep] ];
  945.         pointthickness = Min[0.015, Max[0.024 / factor, 0.003]];
  946.  
  947.         (*  Compute the initial and final zero plot  *)
  948.         If [ ! FreeQ[numer, z],
  949.              {Szplot, Ezplot} =
  950.             initPlot[ (numer /. ktemp -> kmin),
  951.                   (numer /. ktemp -> kmax), z, "O", "*O*" ] ];
  952.  
  953.         (*  Compute the initial and final pole plot  *)
  954.         If [ ! FreeQ[denom, z],
  955.              {Spplot, Epplot} =
  956.             initPlot[ (denom /. ktemp -> kmin),
  957.                   (denom /. ktemp -> kmax), z, "X", "*X*" ] ];
  958.  
  959.         (*  Compute the zero locus  *)
  960.         If [ ! FreeQ[numer, z],
  961.              zplot = oneLocus[ numer, ktemp, kmin, kmax, kstep,
  962.                        z, pointthickness ] ];
  963.  
  964.         (*  Compute the pole locus  *)
  965.         If [ ! FreeQ[denom, z],
  966.              pplot = oneLocus[ denom, ktemp, kmin, kmax, kstep,
  967.                        z, pointthickness ] ];
  968.  
  969.         (*  Enable compilation warnings  *)
  970.         If [ TrueQ[ $VersionNumber >= 2.0 ],
  971.              On[ ParametricPlot::ppcom ] ];
  972.  
  973.         { localplot[ "Pole-Zero Diagram at Start", options,
  974.                  Szplot, Spplot],
  975.           localplot[ "Root Locus", options, Szplot, Spplot,
  976.                  Ezplot, Epplot, zplot, pplot ],
  977.           localplot[ "Pole-Zero Diagram at End", options,
  978.                  Ezplot, Epplot ] }  ]
  979.  
  980.  
  981. (*  D.  P O L E - Z E R O     P L O T T I N G  *)
  982.  
  983. (*      1. S u p p o r t i n g   F u n c t i o n s  *)
  984.  
  985. (*       a.  For graphing the region of convergence as an annulus  *)
  986.  
  987. (*  findline  *)
  988. findline[ycoord_, cos_, sin_] :=
  989.     Block [ {xcoord},
  990.         xcoord = N[ Sqrt[1 - ycoord^2] ];
  991.         { rotate[xcoord, ycoord, cos, sin],
  992.           rotate[-xcoord, ycoord, cos, sin] } ]
  993.  
  994. makeannulus[{rm_, rp_}] := 
  995.     ToCollection[ CirclePS[rp, Dashing[{0.05,0.05}]],
  996.               ShadedAnnulus[rm, rp],
  997.               CirclePS[rm, Dashing[{0.05,0.05}]] ]
  998.  
  999. (*  makefilllines  *)
  1000. (*    Fill lines for a unit circle:  tilt in degrees, *)
  1001. (*  and angular separation ~= 360 / numsamples          *)
  1002. makefilllines[tilt_, numsamples_] :=
  1003.     Block [ {costilt, maxi, sintilt},
  1004.         costilt = N[ Cos[ tilt Pi / 180 ] ];
  1005.         sintilt = N[ Sin[ tilt Pi / 180 ] ];
  1006.         maxi = Floor[numsamples / 4] - 1;
  1007.  
  1008.         Join[ Table[ findline[4 i / numsamples, costilt, sintilt],
  1009.                  {i, 0, maxi} ],
  1010.               Table[ findline[-1 + 4 i / numsamples, costilt, sintilt],
  1011.                  {i, 1, maxi} ] ] ]
  1012.  
  1013. (*  rotate  *)
  1014. rotate[x_, y_, cos_, sin_] := { x cos - y sin, x sin + y cos }
  1015.  
  1016. (*  splitline  *)
  1017. splitline[ap_, rm_, slope_, {point1_, point2_}] :=
  1018.     Block [    {b, sqrtterm, x1, x2, y1, y2},
  1019.         b = Second[point1] - slope First[point1];
  1020.         sqrtterm = Re[ N[ Sqrt[ rm^2 ap - b^2 ] ] ];
  1021.         x1 = ( - slope b + sqrtterm ) / ap;
  1022.         x2 = ( - slope b - sqrtterm ) / ap;
  1023.         y1 = slope x1 + b;
  1024.         y2 = slope x2 + b;
  1025.         {{point1, {x1, y1}}, {{x2, y2}, point2}} ]
  1026.  
  1027. (*  ShadedAnnulus *)
  1028. AngularTilt = 20                  (* default values *)
  1029. NumSamples = 72
  1030. FillLines = makefilllines[AngularTilt, NumSamples]
  1031. FillLineSlope = N[ Tan[AngularTilt Pi / 180] ]      (* slope of all fill lines *)
  1032.  
  1033. ShadedAnnulus[rm_, rp_] :=
  1034.     ShadedAnnulus[ rm, rp, AngularTilt, FillLines, FillLineSlope ]
  1035.  
  1036. ShadedAnnulus[rm_, rp_, tilt_] :=
  1037.     ShadedAnnulus[ rm, rp, tilt, makefilllines[tilt, NumSamples] ]
  1038.  
  1039. ShadedAnnulus[rm_, rp_, tilt_, filllines_] :=
  1040.     ShadedAnnulus[ rm, rp, tilt, filllines, N[Tan[tilt Pi / 180]] ]
  1041.  
  1042. ShadedAnnulus[0, rp_, tilt_, filllines_, slope_] :=
  1043.     Graphics[ Map[Line, rp filllines] ]
  1044.  
  1045. ShadedAnnulus[rm_, rp_, tilt_, filllines_, slope_] :=
  1046.     Block [    {ap, fillpattern, i, intersections,
  1047.          maxi, newlines, numlines, split, xnew},
  1048.  
  1049.         fillpattern = (rp filllines);
  1050.         numlines = Length[filllines];
  1051.  
  1052.         (* Find intersection of lines with rminus circle *)
  1053.         (* Equations:  x^2 + y^2 = r^2 and y = m x + b   *)
  1054.         (* reduces to quadratic a' x^2 + b' x + c', with *)
  1055.         (* a' = (m^2 + 1), b' = 2 m b, c' = b^2 - r^2     *)
  1056.         ap = ( slope^2 + 1 );
  1057.  
  1058.         (* Intersection of first line, zero intercept     *)
  1059.         xnew = N[ - rm / Sqrt[ap] ];
  1060.         newlines = { { {xnew, slope xnew},
  1061.                    Second[ First[fillpattern] ] } };
  1062.         fillpattern[[1]] = { First[ First[fillpattern] ],
  1063.                      {-xnew, - slope xnew} };
  1064.  
  1065.         (* Find intersections of other lines         *)
  1066.         intersections = ( numlines + 2 ) rm / rp;
  1067.         maxi = Floor[ intersections / 2];
  1068.         For [ i = 1, i <= maxi, i++,
  1069.               split = splitline[ ap, rm, slope, fillpattern[[i + 1]] ];
  1070.               AppendTo[ newlines, First[split] ];
  1071.               fillpattern[[i + 1]] = Second[split];
  1072.               split = splitline[ ap, rm, slope,
  1073.                      fillpattern[[numlines - i + 1]] ];
  1074.               AppendTo[ newlines, First[split] ];
  1075.               fillpattern[[numlines - i + 1]] = Second[split] ];
  1076.  
  1077.         Graphics[ Map[Line, Join[fillpattern, newlines]] ] ]
  1078.  
  1079. (*  ROCPlot  *)
  1080. Options[ROCPlot] :=
  1081.   { Axes -> axesdefaultvalue, AspectRatio -> 1,
  1082.     DisplayFunction :> $DisplayFunction }
  1083.  
  1084. ROCPlot[{rm_, Infinity}, options___] :=
  1085.     Block [    {oplist, plotrange, pmax, rmax},
  1086.         oplist = ToList[options];
  1087.         plotrange = Replace[PlotRange, oplist];
  1088.         pmax = If [ SameQ[plotrange, PlotRange],
  1089.                 N[ Sqrt[2] rm ],
  1090.                 Max[ N[Sqrt[2] rm],
  1091.                  Abs[Select[N[Flatten[plotrange]], NumberQ]] ] ];
  1092.         rmax = 1.75 pmax;     (* pmax times some number > Sqrt[2] *)
  1093.         Apply [    Show,
  1094.             Join [ {makeannulus[{rm, rmax}]},
  1095.                    oplist,
  1096.                    Options[ROCPlot],
  1097.                    PlotRange -> {{-pmax, pmax}, {-pmax, pmax}} ] ] ]
  1098.  
  1099. ROCPlot[{rm_, rp_}, options___] :=
  1100.     Apply[ Show,
  1101.            Join[{makeannulus[{rm, rp}]},
  1102.             ToList[options],
  1103.             Options[ROCPlot]] ]
  1104.  
  1105. (*       b.  other supporting functions  *)
  1106.  
  1107.  
  1108.  
  1109. (*  printInvalidPoles  *)
  1110. printInvalidPoles[invalidROClist_, polelist_, rm_, rp_, op_, smallestPole_] :=
  1111.     Block [ {rocs, pairs},
  1112.         Message[PoleZeroPlot::invalidROC, invalidROClist, rm, rp];
  1113.         roclist = Sort[ Map[ op, polelist ] ];
  1114.         pairs = Transpose[ { Prepend[roclist, smallestPole],
  1115.                      Append[roclist, Infinity] } ];
  1116.         Print[ PoleZeroPlot::posROC, Rationalize[pairs] ] ]
  1117.                  
  1118. (*  InvalidZPoleQ  *)
  1119. InvalidZPoleQ[rminus_, r_, rplus_] := TrueQ[rminus < Abs[N[r]] < rplus]
  1120.  
  1121. (*  InvalidSPoleQ  *)
  1122. InvalidSPoleQ[rminus_, r_, rplus_] := TrueQ[rminus < Re[N[r]] < rplus]
  1123.  
  1124. (*  MyRationalPolyQ  *)
  1125. MyRationalPolyQ[f_, z_] :=
  1126.     MixedPolynomialQ[ Expand[Numerator[f]], z ] &&
  1127.     MixedPolynomialQ[ Expand[Denominator[f]], z ]
  1128.  
  1129. (*  pzOptions  *)
  1130. pzOptions[ op___ ] :=
  1131.     Join[ ToList[op],
  1132.           ToList[ Options[PoleZeroPlot] ], 
  1133.           { DisplayFunction :> $DisplayFunction } ]
  1134.  
  1135. (*  SeparableQ  *)
  1136. Separate[x_ y_, {var1_, var2_}] := { x, y } /; FreeQ[x, var2] && FreeQ[y, var1]
  1137. Separable[x_ y_, {var1_, var2_}] := True /; FreeQ[x, var2] && FreeQ[y, var1]
  1138. SeparableQ[p_, vars_] := TrueQ[ Separable[p, vars] ]
  1139.  
  1140. (*      2.  D r i v e r  *)
  1141.  
  1142. Options[ PoleZeroPlot ] :=
  1143.     { Dialogue -> True, DisplayFunction :> $DisplayFunction }
  1144.  
  1145. (*  One-dimensional pole-zero plotting for z-transform objects  *)
  1146. PoleZeroPlot[ ZTransData[f_, Rminus[rm_], Rplus[rp_], ZVariables[z_]], op___ ] :=
  1147.     If [ Length[z] <= 2,
  1148.          PoleZeroPlot[ f, z, rm, rp, True, op ],
  1149.          Message[PoleZeroPlot::multidimensions] ]
  1150.  
  1151. (*  One-dimensional pole-zero plotting for Laplace transform objects  *)
  1152. PoleZeroPlot[ LTransData[f_, Rminus[rm_], Rplus[rp_], LVariables[s_]], op___ ] :=
  1153.     If [ Length[s] <= 2,
  1154.          PoleZeroPlot[ f, s, rm, rp, False, op ],
  1155.          Message[PoleZeroPlot::multidimensions] ]
  1156.  
  1157. (*  Two-dimensional pole-zero plotting driver                 *)
  1158. (*  --  separable transforms generate two separable pole-zero plots  *) 
  1159. (*  --  symmetric transfer function f only requires one plot         *)
  1160. (*  --  otherwise, project z1 onto z2 and z2 onto z1             *)
  1161. PoleZeroPlot[f_, {z1_Symbol, z2_Symbol}, {rm1_, rm2_}, {rp1_, rp2_}, rest__ ] :=
  1162.     Block [ {},
  1163.         Needs[ "SignalProcessing`Support`MDPlotting`" ];
  1164.         poleZeroPlot2D[f, {z1, z2}, {rm1, rm2}, {rp1, rp2}, rest ] ]
  1165.  
  1166.  
  1167. (*      3.  P l o t t i n g   R o u t i n e s   f o r   z - d o m a i n  *)
  1168.  
  1169. PoleZeroPlot[f_, z_Symbol, rm_, rp_, True, op___] :=    
  1170.     Block [ {abspolelist, fillplot, invalidROClist, options,
  1171.          polelist, poleplot, rmax, rminus, rplus, rminuscircle,
  1172.          rpluscircle, unitcircle, xlabel, ylabel, zeroplot,
  1173.          zstr, ztrans},
  1174.  
  1175.         options = pzOptions[op];
  1176.  
  1177.         rmax = 1;                  (* default values *)
  1178.         unitcircle = CirclePS[1];
  1179.  
  1180.         rminus = N[rm];              (* numeric approximation *)
  1181.         rplus = N[rp];
  1182.  
  1183.         (*  Analyze the z-transform function  *)
  1184.  
  1185.         ztrans = Together[f];
  1186.         If [ ! MyRationalPolyQ[ztrans, z],
  1187.              Message[ PoleZeroPlot::notrational ];
  1188.              Message[ PoleZeroPlot::noplot ];
  1189.              Return[{}] ];
  1190.  
  1191.         (*  Find the poles and zeroes  *)
  1192.  
  1193.         zerolist = GetRootList[Numerator[ztrans], z];
  1194.         polelist = GetRootList[Denominator[ztrans], z];
  1195.         If [ InformUserQ[options],
  1196.              Print[ " " ];
  1197.              Print[ "The zeroes are:  ", zerolist ];
  1198.              Print[ "The poles are:   ", polelist ] ];
  1199.  
  1200.         (*  Check consistency of poles with region of convergence  *)
  1201.         (*  Exit function if inconsistency encountered           *)
  1202.  
  1203.         invalidROClist = Select[ polelist,
  1204.                      InvalidZPoleQ[rminus, #1, rplus]& ];
  1205.         If [ ! SameQ[invalidROClist, {}],
  1206.              printInvalidPoles[ invalidROClist, polelist,
  1207.                     rm, rp, Abs, 0 ];
  1208.              Return[polelist] ];
  1209.  
  1210.         (*  Compute plot of poles and zeroes        *)
  1211.         (*  Determine the plot's its circular extent    *)
  1212.  
  1213.         rmax = Max[ Join[ {rmax},
  1214.                   Map[Abs, ToList[N[zerolist]]],
  1215.                   Map[Abs, ToList[N[polelist]]] ] ];
  1216.  
  1217.         zeroplot = If [ EmptyQ[zerolist],
  1218.                          NullPlot,
  1219.                 PointwisePlot[ComplexTo2DCoordList[zerolist],
  1220.                           "O", "*O*"] ];
  1221.  
  1222.         poleplot = If [ EmptyQ[polelist],
  1223.                 NullPlot,
  1224.                 PointwisePlot[ComplexTo2DCoordList[polelist],
  1225.                           "X", "*X*"] ];
  1226.  
  1227.         (*  Build region of convergence plot  *)
  1228.  
  1229.         If [ NumberQ[rminus],
  1230.              rminuscircle = CirclePS[ rminus, Dashing[{0.05,0.05}] ];
  1231.              rmax = Max[rmax, Abs[rminus]],
  1232.              rminuscircle = NullPlot ];
  1233.  
  1234.         If [ NumberQ[rplus],
  1235.              rpluscircle = CirclePS[ rplus, Dashing[{0.05,0.05}] ];
  1236.              rmax = Max[rmax, Abs[rplus]],
  1237.              rpluscircle = NullPlot ];
  1238.  
  1239.         rplus = N[ If[ SameQ[rp, Infinity], 2 rmax, rp] ];  (* kludge *)
  1240.  
  1241.         fillplot = If [ NumberQ[rminus] && NumberQ[rplus],
  1242.                 ShadedAnnulus[rminus, rplus],
  1243.                 NullPlot ];
  1244.  
  1245.         (*  Send the contour plot to the screen  *)
  1246.         rmax *= 1.25;
  1247.         zstr = StripPackage[z];
  1248.         xlabel = "Re " ~StringJoin~ zstr;
  1249.         ylabel = "Im " ~StringJoin~ zstr;
  1250.  
  1251.         options = Join[ RemoveOptions[options, {Dialogue}];
  1252.                 { PlotRange -> {{-rmax, rmax}, {-rmax, rmax}},
  1253.                   AspectRatio -> 1,
  1254.                   AxesLabel -> { xlabel, ylabel },
  1255.                   Axes -> axesdefaultvalue } ];
  1256.  
  1257.         Show [unitcircle, rminuscircle, rpluscircle,
  1258.               fillplot, zeroplot, poleplot, options];
  1259.  
  1260.         polelist ]
  1261.  
  1262.  
  1263. (*      4.  P l o t t i n g   R o u t i n e s   f o r   s - d o m a i n  *)
  1264.  
  1265. (*  One-dimensional pole-zero plotting  *)
  1266. PoleZeroPlot[f_, s_Symbol, rm_, rp_, False, op___] :=
  1267.     Block [ {invalidROClist, jj, ltrans, numericlist, options,
  1268.          points, rminus, rminusline, rplus, rplusline, polelist,
  1269.          poleplot, shading, skew, sstr, x1, x2, xlabel, xmin,
  1270.          xmax, w, ylabel, ymin, ymax, ystep, zerolist, zeroplot},
  1271.  
  1272.         options = pzOptions[op];
  1273.  
  1274.         ltrans = Together[f];            (* default values *)
  1275.         rminusline = rplusline = shading = NullPlot;
  1276.  
  1277.         rminus = N[rm];             (* numeric approximation *)
  1278.         rplus = N[rp];
  1279.  
  1280.         If [ ! RationalPolynomialQ[ltrans, s],
  1281.              Message[ PoleZeroPlot::notrational ];
  1282.              Message[ PoleZeroPlot::noplot ];
  1283.              Return[{}] ];
  1284.  
  1285.         zerolist = GetRootList[Numerator[ltrans], s];
  1286.         polelist = GetRootList[Denominator[ltrans], s];
  1287.         If [ InformUserQ[options],
  1288.              Print[ " " ];
  1289.              Print[ "The zeroes are:  ", zerolist ];
  1290.              Print[ "The poles are:   ", polelist ] ];
  1291.  
  1292.         (*  Check consistency of poles with region of convergence  *)
  1293.         (*  Exit function if inconsistency encountered           *)
  1294.  
  1295.         invalidROClist = Select[ polelist,
  1296.                      InvalidSPoleQ[rminus, #1, rplus]& ];
  1297.         If [ ! SameQ[invalidROClist, {}],
  1298.              printInvalidPoles[ invalidROClist, polelist,
  1299.                     rm, rp, Re, -Infinity ];
  1300.              Return[polelist] ];
  1301.  
  1302.         numericlist = ToList[N[zerolist]] ~Join~ ToList[N[polelist]];
  1303.  
  1304.         (*  Determine extent of pole-zero plot        *)
  1305.         (*  Mathematica does not do this automatically    *)
  1306.  
  1307.         xmin = Min[Re[numericlist]];
  1308.         xmax = Max[Re[numericlist]];
  1309.  
  1310.         If [ NumberQ[rminus], xmin = Min[xmin, rminus] ];
  1311.         If [ NumberQ[rplus],  xmax = Max[xmax, rplus]  ];
  1312.  
  1313.         xmin = ( 1.0 - 0.25 Sign[xmin] ) xmin;    (* move toward -Inf *)
  1314.         xmax = ( 1.0 + 0.25 Sign[xmax] ) xmax;    (* move toward +Inf *)
  1315.         If [ ZeroQ[xmax], xmax = -.1 xmin ];
  1316.  
  1317.         ymin = Min[Im[numericlist]];
  1318.         ymax = Max[Im[numericlist]];
  1319.  
  1320.         If [ SameQ[ymin, ymax],
  1321.              ymin = ymin - 1;
  1322.                ymax = ymax + 1,
  1323.              ymin = ( 1.0 - 0.25 Sign[ymin] ) ymin;
  1324.                ymax = ( 1.0 + 0.25 Sign[ymax] ) ymax ];
  1325.  
  1326.         (*  Build convergence boundaries (lines)  *)
  1327.  
  1328.         If [ NumberQ[rminus],
  1329.              rminusline = Graphics[ { Dashing[{0.05, 0.05}],
  1330.                           Line[{{rminus, ymin},
  1331.                             {rminus, ymax}}] } ] ];
  1332.  
  1333.         If [ NumberQ[rplus],
  1334.              rplusline = Graphics[ { Dashing[{0.05, 0.05}],
  1335.                          Line[{{rplus, ymin},
  1336.                            {rplus, ymax}}] } ] ];
  1337.  
  1338.         If [ (NumberQ[rminus] || SameQ[rm, -Infinity]) &&
  1339.                (NumberQ[rplus] || SameQ[rp, Infinity]),
  1340.              x1 = Max[rminus, xmin];
  1341.              x2 = Min[rplus, xmax];
  1342.              ystep = ( ymax - ymin ) / 15;
  1343.              skew = 4 ystep;
  1344.              points = Table[ Line[{{x1,jj}, {x2,jj+skew}}],
  1345.                      {jj, ymin - skew, ymax, ystep} ];
  1346.              shading = Graphics[ points ] ];
  1347.  
  1348.         (*  Build plots of zeroes as O's and poles as X's  *)
  1349.  
  1350.         If [ EmptyQ[zerolist],
  1351.              zeroplot = NullPlot,
  1352.              zeroplot = PointwisePlot[ ComplexTo2DCoordList[zerolist],
  1353.                            "O", "*O*" ] ];
  1354.  
  1355.         If [ EmptyQ[polelist],
  1356.              poleplot = NullPlot,
  1357.              poleplot = PointwisePlot[ ComplexTo2DCoordList[polelist],
  1358.                            "X", "*X*" ] ];
  1359.  
  1360.         (*  Plot pole-zero diagram, all together now  *)
  1361.         sstr = StripPackage[s];
  1362.         xlabel = "Re " ~StringJoin~ sstr;
  1363.         ylabel = "Im " ~StringJoin~ sstr;
  1364.  
  1365.         options = Join[ RemoveOptions[options, {Dialogue}];
  1366.                 { PlotRange -> {{xmin, xmax}, {ymin, ymax}},
  1367.                   AxesLabel -> { xlabel, ylabel },
  1368.                   Axes -> axesdefaultvalue } ];
  1369.  
  1370.         Show [ rminusline, rplusline, shading,
  1371.                zeroplot, poleplot, options ];
  1372.  
  1373.         polelist ]
  1374.  
  1375.  
  1376.  
  1377. (*  E.  S T A B I L I T Y     A N A L Y S I S  *)
  1378.  
  1379. (*
  1380.       The trick here for m-D transforms is to replace any occurrence of
  1381.     the transform variables in the ROC with the value of the boundary at
  1382.     which the ROC becomes unstable.  For the Laplace transform, this is
  1383.     the line Re(s) = 0.  For the z-transform, this is the m-D unit circle;
  1384.     because the ROC in the z-transform are real numbers (absolute values
  1385.     of complex-valued quantities), we let z1 = 1, z2 = 1, etc., instead
  1386.     of z1 = Exp[I w1], z2 = Exp[I w2], etc.
  1387.  *)
  1388.  
  1389. unknown[ condexpr_ ] :=
  1390.         Block [ {simplified},
  1391.                 simplified = Simplify[condexpr] /. SPLessGreaterRules;
  1392.                 If [ N[simplified], True, False, simplified ] ]
  1393.  
  1394. LTransData/: Stable[ LTransData[t_, rm_, rp_, v_] ] :=
  1395.         Block [ {inrange, vars, ltrans},
  1396.                 ltrans = LTransData[t, rm, rp, v];
  1397.                 inrange = InRange[GetRMinus[ltrans], 0, GetRPlus[ltrans],
  1398.                                   -Infinity, Infinity, Less, Less];
  1399.                 vars = First[v];
  1400.                 If [ ListQ[vars],           (* evaluate along imag. axes *)
  1401.                      inrange = inrange /. ( Map[Re, vars] ~ReplaceWith~ 0 ) ];
  1402.                 If [ inrange, True, False, unknown[inrange] ] ]
  1403.  
  1404. ZTransData/: Stable[ ZTransData[t_, rm_, rp_, v_] ] :=
  1405.         Block [ {inrange, vars, ztrans},
  1406.                 ztrans = ZTransData[t, rm, rp, v];
  1407.                 inrange = InRange[GetRMinus[ztrans], 1, GetRPlus[ztrans],
  1408.                                   0, Infinity, Less, Less];
  1409.                 vars = First[v];
  1410.                 If [ ListQ[vars],           (* evaluate along unit circles *)
  1411.                      inrange = inrange /. ( vars ~ReplaceWith~ 1 ) ];
  1412.                 If [ inrange, True, False, unknown[inrange] ] ]
  1413.  
  1414.  
  1415.  
  1416.  
  1417. (*  E N D     P A C K A G E  *)
  1418.  
  1419. End[]
  1420. EndPackage[]
  1421.  
  1422.  
  1423. (*  A L I A S E S  *)
  1424.  
  1425. Unprotect[ MagnitudePhasePlot ]
  1426. Clear[ MagnitudePhasePlot ]
  1427. MagnitudePhasePlot = MagPhasePlot
  1428. MagnitudePhasePlot::usage = MagPhasePlot::usage
  1429. Protect[ MagnitudePhasePlot ]
  1430.  
  1431. Unprotect[ RootLocusPlot ]
  1432. Clear[ RootLocusPlot ]
  1433. RootLocusPlot = RootLocus
  1434. RootLocusPlot::usage = RootLocus::usage
  1435. Protect[ RootLocusPlot ]
  1436.  
  1437. If [ TrueQ[ $VersionNumber >= 2.0 ],
  1438.      On[ General::spell ];
  1439.      On[ General::spell1 ] ]
  1440.  
  1441.  
  1442. (*  H E L P     I N F O R M A T I O N  *)
  1443.  
  1444. Block [    {newfuns},
  1445.     newfuns =
  1446.       { AddT,        ConjT,            ConvolveT,
  1447.         DerivativeT,    IntegrateT,        LineImpulsemDT,
  1448.         MagPhasePlot,    MultiDInvTransform,    MultiDLookup,
  1449.         MultiDTransform,    MultT,            PoleZeroPlot,
  1450.         ROCPlot,        RootLocus,        ScaleT,
  1451.         ShadedAnnulus,    Stable,            SubstituteForT,
  1452.         SubT,        TransformDialogue,    TransformFixUp };
  1453.     Combine[ SPfunctions, newfuns ];
  1454.     Apply[ Protect, newfuns ] ]
  1455.  
  1456.  
  1457. (*  E N D I N G     M E S S A G E  *)
  1458.  
  1459. Print["Supporting routines for transform rule bases are loaded."]
  1460. Null
  1461.